# Clear
rm(list = ls())

# Working Directory
getwd()
list.files()
setwd("/Users/claywebb/Dropbox/KU/Research/Sanctions/Power Politics or Public Pandering/IIRR")

## Load Dataset
sanctions <- ts(read.csv("sanc_rep_data.csv"), start = c(1978,1), end = c(2005,12), frequency = 12) 

###################################################
### Section 2.1: Contemporaneous Identification ###
###################################################

# Load MSBVAR Package
library(MSBVAR)

# Time Series
VAR_SER <- sanctions[,c("sanc","approval","unem1","dt_inflation","umics")]

# Deterministic Variables
z1 <- ts(as.matrix(cbind(sanctions[,"elecyear"],
                         sanctions[,"inauguration"],
                         sanctions[,"reagan"],
                         sanctions[,"bush1"],
                         sanctions[,"clinton"],
                         sanctions[,"bush2"],
                         sanctions[,"iranhostage"],
                         sanctions[,"reaganshot"],
                         sanctions[,"grenada"],
                         sanctions[,"libyabomb"],
                         sanctions[,"irancontra"],
                         sanctions[,"panamainvade"],
                         sanctions[,"budgetsummit"],
                         sanctions[,"persiangulf"],
                         sanctions[,"somalia"],
                         sanctions[,"haitiinvade"],
                         sanctions[,"lewinsky"],
                         sanctions[,"kosovoinvade"],
                         sanctions[,"iraqifreedom"],
                         sanctions[,"chinaspyplane"],
                         sanctions[,"afghanwar"],
                         sanctions[,"sep11"], 
                         sanctions[,"iraqinvade"], 
                         sanctions[,"unscom"])), start = c(1978,1), end = c(2005,12), frequency = 12) 
z1

## Set Lag Length 
p <- 2

## Set Hyperparameters
lambda0 <- 0.5
lambda1 <- 0.1
lambda3 <- 2
lambda4 <- 0.5
lambda5 <- 0.05
mu5 <- 1
mu6 <- 0
qm <- 12

m <- ncol(VAR_SER)

## 1. A0 Test
test.A0 <- diag(m)
colnames(test.A0) <- colnames(VAR_SER)
rownames(test.A0) <- colnames(VAR_SER)
test.A0
test.A0[2,1] <- 1
test.A0[c(1),2] <- 1
test.A0

## 2. A0 Sanctions
sanctions.A0 <- diag(m)
colnames(sanctions.A0) <- colnames(VAR_SER)
rownames(sanctions.A0) <- colnames(VAR_SER)
sanctions.A0
sanctions.A0[2,1] <- 1
sanctions.A0[c(1),2] <- 1
sanctions.A0[c(1),3] <- 1
sanctions.A0[c(1),4] <- 1
sanctions.A0[c(1),5] <- 1
sanctions.A0

## 3. A0 Approval
approval.A0 <- diag(m)
colnames(approval.A0) <- colnames(VAR_SER)
rownames(approval.A0) <- colnames(VAR_SER)
approval.A0
approval.A0[2,1] <- 1
approval.A0[c(1),2] <- 1
approval.A0[c(2),3] <- 1
approval.A0[c(2),4] <- 1
approval.A0[c(2),5] <- 1
approval.A0

## 4. A0 Sanctions Cause and Effect
sanctions_ce.A0 <- diag(m)
colnames(sanctions_ce.A0) <- colnames(VAR_SER)
rownames(sanctions_ce.A0) <- colnames(VAR_SER)
sanctions_ce.A0
sanctions_ce.A0[c(2,3,4,5),1] <- 1
sanctions_ce.A0[c(1),2] <- 1
sanctions_ce.A0[c(1),3] <- 1
sanctions_ce.A0[c(1),4] <- 1
sanctions_ce.A0[c(1),5] <- 1
sanctions_ce.A0

## 5. A0 Approval Cause and Effect
approval_ce.A0 <- diag(m)
approval_ce.A0
approval_ce.A0[c(2),1] <- 1
approval_ce.A0[c(1,3,4,5),2] <- 1
approval_ce.A0[c(2),3] <- 1
approval_ce.A0[c(2),4] <- 1
approval_ce.A0[c(2),5] <- 1
colnames(approval_ce.A0) <- colnames(VAR_SER)
rownames(approval_ce.A0) <- colnames(VAR_SER)
approval_ce.A0

## 6. A0_sanc_app_e
sanc_app_e.A0 <- diag(m)
sanc_app_e.A0
sanc_app_e.A0[c(2),1] <- 1
sanc_app_e.A0[c(1),2] <- 1
sanc_app_e.A0[c(1,2,5),3] <- 1
sanc_app_e.A0[c(1,2,5),4] <- 1
sanc_app_e.A0[c(1,2),5] <- 1
sanc_app_e.A0
colnames(sanc_app_e.A0) <- colnames(VAR_SER)
rownames(sanc_app_e.A0) <- colnames(VAR_SER)
sanc_app_e.A0

## 7. A0_sanc_app_c
sanc_app_c.A0 <- diag(m)
sanc_app_c.A0
sanc_app_c.A0[c(2,3,4,5),1] <- 1
sanc_app_c.A0[c(1,3,4,5),2] <- 1
sanc_app_c.A0
colnames(sanc_app_c.A0) <- colnames(VAR_SER)
rownames(sanc_app_c.A0) <- colnames(VAR_SER)
sanc_app_c.A0

## 8. Econ_A0
econ.A0 <- diag(m)
econ.A0
econ.A0[c(2,3,4,5),1] <- 1
econ.A0[c(1,5),2] <- 1
econ.A0[c(4,5),3] <- 1
econ.A0[c(3,5),4] <- 1
colnames(econ.A0) <- colnames(VAR_SER)
rownames(econ.A0) <- colnames(VAR_SER)
econ.A0

## 9. Pol
pol.A0 <- diag(m)
colnames(pol.A0) <- colnames(VAR_SER)
rownames(pol.A0) <- colnames(VAR_SER)
pol.A0
pol.A0[c(2,3,4,5),1] <- 1
pol.A0[c(1,5),2] <- 1
pol.A0[c(5),3] <- 1
pol.A0[c(5),4] <- 1
pol.A0[c(2),5] <- 1
pol.A0

## 10. A0_Recursive
ident <- matrix(0, m, m)
for(i in 0:(m-1)) ident[,i] <- c(rep(0, i), rep(1,m-i))
diag(ident) <- 1
ident
recursive.A0 <- ident
colnames(recursive.A0) <- colnames(VAR_SER)
rownames(recursive.A0) <- colnames(VAR_SER)
recursive.A0

## Estimate Models
test_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                    lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                    mu5 = mu5, mu6 = mu6, qm = qm, test.A0)

sanctions_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                         lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                         mu5 = mu5, mu6 = mu6, qm = qm, sanctions.A0)

approval_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                        lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                        mu5 = mu5, mu6 = mu6, qm = qm, approval.A0)

sanctions_ce_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                            lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                            mu5 = mu5, mu6 = mu6, qm = qm, sanctions_ce.A0)

approval_ce_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                           lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                           mu5 = mu5, mu6 = mu6, qm = qm, approval_ce.A0)

sanc_app_e_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                          lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                          mu5 = mu5, mu6 = mu6, qm = qm, sanc_app_e.A0)

sanc_app_c_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                          lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                          mu5 = mu5, mu6 = mu6, qm = qm, sanc_app_c.A0)

econ_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                    lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                    mu5 = mu5, mu6 = mu6, qm = qm, econ.A0)

pol_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                   lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                   mu5 = mu5, mu6 = mu6, qm = qm, pol.A0)

recursive_mod <- szbsvar(Y = VAR_SER, p = p, z = z1, 
                         lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                         mu5 = mu5, mu6 = mu6, qm = qm, recursive.A0)

# Sample

## Burnin 
N1 <- 10000

## Gibbs Draws
N2 <- 30000

## Thin Integer
thin <- 2

## Test Mod

test_mod_out <- gibbs.A0(varobj = test_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

sanctions_mod_out <- gibbs.A0(varobj = sanctions_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

approval_mod_out <- gibbs.A0(varobj = approval_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

sanctions_ce_mod_out <- gibbs.A0(varobj = sanctions_ce_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

approval_ce_mod_out <- gibbs.A0(varobj = approval_ce_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

sanc_app_e_mod_out <- gibbs.A0(varobj = sanc_app_e_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

sanc_app_c_mod_out <- gibbs.A0(varobj = sanc_app_c_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

econ_mod_out <- gibbs.A0(varobj = econ_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

pol_mod_out <- gibbs.A0(varobj = pol_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

recursive_mod_out <- gibbs.A0(varobj = recursive_mod, N1 = N1, N2 = N2, thin = thin, normalization = "DistanceMLA")

test_mod_post <- posterior.fit(test_mod, test_mod_out, maxiterbs = 500)
sanctions_mod_post <- posterior.fit(sanctions_mod, sanctions_mod_out, maxiterbs = 500)
approval_mod_post <- posterior.fit(approval_mod, approval_mod_out, maxiterbs = 500)
sanctions_ce_mod_post <- posterior.fit(sanctions_ce_mod, sanctions_ce_mod_out, maxiterbs = 500)
approval_ce_mod_post <- posterior.fit(approval_ce_mod, approval_ce_mod_out, maxiterbs = 500)
sanc_app_e_mod_post <- posterior.fit(sanc_app_e_mod, sanc_app_e_mod_out, maxiterbs = 500)
sanc_app_c_mod_post <- posterior.fit(sanc_app_c_mod, sanc_app_c_mod_out, maxiterbs = 500)
econ_mod_post <- posterior.fit(econ_mod, econ_mod_out, maxiterbs = 500)
pol_mod_post <- posterior.fit(pol_mod, pol_mod_out, maxiterbs = 500)
recursive_mod_post <- posterior.fit(recursive_mod, recursive_mod_out, maxiterbs = 500)

## Function to Organize Posterior Fit Output
pf_org <- function(fits_ob){
  outvec <- c(fits_ob$log.prior, fits_ob$log.llf, fits_ob$log.posterior.Aplus, sum(fits_ob$log.marginal.A0k), fits_ob$log.marginal.data.density)
  outvec
}

## Table 3
tab_3 <- round(rbind(pf_org(test_mod_post),
pf_org(sanctions_mod_post),
pf_org(approval_mod_post),
pf_org(sanctions_ce_mod_post),
pf_org(approval_ce_mod_post),
pf_org(sanc_app_c_mod_post),
pf_org(sanc_app_e_mod_post),
pf_org(econ_mod_post),
pf_org(pol_mod_post),
pf_org(recursive_mod_post)),3)

rownames(tab_3) <- c("Test Model","Sanctions","Approval","Sanctions Cause and Effect","Approval Cause and Effect", 
                     "Sanctions and Approvl Cause", "Sanctions and Approval Effect", "Economics", "Politics","Recursive")
tab_3


#########################################
### Section 1.3: Lag Length Selection ###
#########################################
lag_sel <- var.lag.specification(VAR_SER,5)
round(lag_sel$ldets,3)
round(lag_sel$results,3)

############################
### Section 2: Variables ###
############################
tab_sanctions <- rbind(round(c(length(sanctions[,"sanc"]),mean(sanctions[,"sanc"]),sd(sanctions[,"sanc"]),min(sanctions[,"sanc"]),max(sanctions[,"sanc"])),2),
      round(c(length(sanctions[,"approval"]),mean(sanctions[,"approval"]),sd(sanctions[,"approval"]),min(sanctions[,"approval"]),max(sanctions[,"approval"])),2),
      round(c(length(sanctions[,"unem1"]),mean(sanctions[,"unem1"]),sd(sanctions[,"unem1"]),min(sanctions[,"unem1"]),max(sanctions[,"unem1"])),2),
      round(c(length(sanctions[,"dt_inflation"]),mean(sanctions[,"dt_inflation"]),sd(sanctions[,"dt_inflation"]),min(sanctions[,"dt_inflation"]),max(sanctions[,"dt_inflation"])),2),
      round(c(length(sanctions[,"umics"]),mean(sanctions[,"umics"]),sd(sanctions[,"umics"]),min(sanctions[,"umics"]),max(sanctions[,"umics"])),2))

rownames(tab_sanctions) <- c("Sanctions","Approval","Unemployment","Inflation","Consumer Sentiment")
colnames(tab_sanctions) <- c("Observations","Mean","SD","Min","Max")
tab_sanctions

#########################################
### Section 3: Posterior Coefficients ###
#########################################
set.seed(121212)

## Estimate Sructural VAR
sanc_app_e_mod <- szbsvar(Y = VAR_SER, z = z1,
                          p = p, lambda0 = lambda0, lambda1 = lambda1, lambda3 = lambda3,  lambda4 = lambda4, lambda5 = lambda5,  
                          mu5 = mu5, mu6 = mu6, qm = qm, ident = sanc_app_e.A0)


pos_coefs <- round(sanc_app_e_mod$B.posterior,2)
colnames(pos_coefs) <- c("Sanctions","Approval","Unemployment","Inflation","ICS")
rownames(pos_coefs) <- c("$\text{Sanctions}_{t-1}$","$\text{Sanctions}_{t-2}$","$\text{Approval}_{t-1}$","$\text{Approval}_{t-2}$",
                         "$\text{Unemployment}_{t-1}$","$\text{Unemployment}_{t-2}$","$\text{Inflation}_{t-1}$","$\text{Inflation}_{t-2}$",
                         "$\text{UMICS}_{t-1}$","$\text{UMICS}_{t-2}$","Constant","Election Year","Inaugruation", "Reagan","H.W. Bush","Clinton",
                         "W Bush", "Iran Hostage", "Reagan Shot", "Grenada", "Libya Bombing", "Iran Contra", "Panama Invasion", "Budget Summit",
                         "Persian Gulf", "Somalia", "Haiti", "Lewinksy Scandal", "Kosovo", "Iraqi Freedom", "China Spy Plane", "Afghan War", "9/11",
                         "Iraq Invasions", "Iraq Bombing")
pos_coefs


##############################
### Section 4: Issue Plots ###
##############################
pol_inf.full <- dynlm::dynlm(approval ~ L(approval, 1) + pol_inf +  unem1 + inflation + umics + 
                               elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                               reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                               persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                               chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

con_mil.full <- dynlm::dynlm(approval ~ L(approval, 1) + con_mil +  unem1 + inflation + umics + 
                               elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                               reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                               persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                               chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

destab_reg.full <- dynlm::dynlm(approval ~ L(approval, 1) + destab_reg +  unem1 + inflation + umics + 
                                  elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                                  reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                                  persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                                  chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

rel_cpm.full <- dynlm::dynlm(approval ~ L(approval, 1) + rel_cpm +  unem1 + inflation + umics + 
                               elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                               reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                               persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                               chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

territory.full <- dynlm::dynlm(approval ~ L(approval, 1) + territory +  unem1 + inflation + umics + 
                                 elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                                 reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                                 persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                                 chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

strat_mat.full <- dynlm::dynlm(approval ~ L(approval, 1) + strat_mat +  unem1 + inflation + umics + 
                                 elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                                 reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                                 persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                                 chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

ally.full <- dynlm::dynlm(approval ~ L(approval, 1) + ally +  unem1 + inflation + umics + 
                            elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                            reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                            persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                            chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

hr.full <- dynlm::dynlm(approval ~ L(approval, 1) + hr +  unem1 + inflation + umics + 
                          elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                          reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                          persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                          chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

prolif.full <- dynlm::dynlm(approval ~ L(approval, 1) + prolif +  unem1 + inflation + umics + 
                              elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                              reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                              persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                              chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

nonstate.full <- dynlm::dynlm(approval ~ L(approval, 1) + nonstate +  unem1 + inflation + umics + 
                                elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                                reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                                persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                                chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

drug.full <- dynlm::dynlm(approval ~ L(approval, 1) + drug +  unem1 + inflation + umics + 
                            elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                            reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                            persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                            chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

environment.full <- dynlm::dynlm(approval ~ L(approval, 1) + environment +  unem1 + inflation + umics + 
                                   elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                                   reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                                   persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                                   chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

trade.full <- dynlm::dynlm(approval ~ L(approval, 1) + trade +  unem1 + inflation + umics + 
                             elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                             reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                             persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                             chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

econ.full <- dynlm::dynlm(approval ~ L(approval, 1) + econ +  unem1 + inflation + umics + 
                            elecyear + inauguration + reagan + bush1 + clinton + bush2 + iranhostage +
                            reaganshot + grenada + libyabomb + irancontra + panamainvade + budgetsummit +
                            persiangulf + somalia + haitiinvade + lewinsky + kosovoinvade + iraqifreedom +
                            chinaspyplane + afghanwar + sep11 + iraqinvade + unscom, data = sanctions)

library(coefplot)
#devtools::install_github("jaredlander/coefplot")
issue_plot <- coefplot::multiplot(econ.full,
                                  trade.full,
                                  environment.full,
                                  drug.full,
                                  nonstate.full,
                                  prolif.full,
                                  hr.full, 
                                  ally.full,
                                  strat_mat.full,
                                  territory.full,
                                  rel_cpm.full,
                                  destab_reg.full,
                                  con_mil.full,
                                  pol_inf.full,
                                  coefficients = c("econ",
                                                   "trade",
                                                   "environment",
                                                   "drug",
                                                   "nonstate",
                                                   "prolif",
                                                   "hr",
                                                   "ally",
                                                   "strat_mat",
                                                   "territory",
                                                   "rel_cpm",
                                                   "destab_reg",
                                                   "con_mil",
                                                   "pol_inf"),
                                  title = "",
                                  xlab = "",
                                  ylab = "",
                                  innerCI = 2,
                                  lwdInner = .25,
                                  zeroColor = "grey",
                                  zeroLWD = 2,
                                  by="Coefficient",
                                  newNames = c(econ = "Implement Economic Reform",
                                               trade = "Trade Practices",
                                               environment = "Improve Environmental Policies",
                                               drug = "Deter or Punish Drug Trafficking Practices",
                                               nonstate = "Terminate Support of Non-State Actors",
                                               prolif = "End Weapons / Material Proliferation",
                                               hr = "Improve Human Rights",
                                               ally = "Retaliate for Alliance or Alignment Choice",
                                               strat_mat = "Deny Strategic Materials",
                                               territory = "Solve Territorial Dispute",
                                               rel_cpm = "Release Citizens, Property, or Material",
                                               destab_reg = "Destabilize Regime",
                                               con_mil = "Contain Military Behavior",
                                               pol_inf = "Contain Political Influence")) + 
  theme(legend.position = "none") + 
  theme(panel.background = element_rect(fill = "white")) +
  theme(axis.text = element_text(colour = "black")) +
  theme(axis.ticks = element_line(color = "black")) +
  geom_vline(xintercept = -15) +
  geom_vline(xintercept = -10) +
  geom_vline(xintercept = -5) +
  geom_vline(xintercept = 5) +
  geom_vline(xintercept = 10) +
  geom_vline(xintercept = 15) + 
  scale_colour_manual(values=c("black", 
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black",
                               "black")) 
issue_plot
